home *** CD-ROM | disk | FTP | other *** search
/ MacFormat 1994 August / August CD.bin / Shareware / Education / numericalmethods Folder / chap_2 / a2_9.m < prev    next >
Encoding:
Text File  |  1994-06-05  |  4.1 KB  |  175 lines  |  [MATS/MATL]
  1. echo off;
  2. % NUMERICAL METHODS: MATLAB Programs, (c) John H. Mathews 1994
  3. % To accompany the text:
  4. % NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992
  5. % Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A.
  6. % This free software is complements of the author.
  7.  
  8. % Algorithm 2.9 (Nonlinear Seidel Iteration).
  9. % Section    2.6,    Iteration for Nonlinear Systems, Page 108
  10. echo on; clc; format long;  hold off; clear
  11.  
  12. % This program implements fixed point iteration
  13.  
  14. % in 2-dimensions for solving  X = G(X).
  15.  
  16. %    Use the vector notation X = (x y).
  17.  
  18. %    Define and store the functions g1(X) and g2(X) and
  19.  
  20. %    define and store the function   G(X) in the
  21.  
  22. %    M-files   g1.m   g2.m   and   G.m   respectively.
  23.  
  24. pause    % Press any key to continue.
  25.  
  26. clc;
  27. % function z = g1(x,y)
  28. % z = (x.^2 - y + 0.5)./2;
  29.  
  30. % function z = g2(x,y)
  31. % z = (- x.^2 - 4.*y.^2 + 8.*y + 4)./8;
  32.  
  33. % function Z = G(X)
  34. % x = X(1); y = X(2);
  35. % Z = [g1(x,y)  g2(x,y)];
  36.  
  37. delete g1.m
  38. diary g1.m; disp('function z = g1(x,y)');...
  39.             disp('z = (x.^2 - y + 0.5)./2;');...
  40. diary off;
  41.  
  42. delete g2.m
  43. diary g2.m; disp('function z = g2(x,y)');...
  44.             disp('z = (- x.^2 - 4.*y.^2 + 8.*y + 4)./8;');...
  45. diary off;
  46.  
  47. delete G.m
  48. diary G.m; disp('function Z = G(X)');...
  49.            disp('x = X(1); y = X(2);');...
  50.            disp('Z = [g1(x,y)  g2(x,y)];');...
  51. diary off;
  52.  
  53. % Remark. g1.m, g2.m, G.m and fix2dim.m are used for Algorithm 2.9
  54. g1(0,0); g2(0,0); G([0 0]); % Test for files g1.m g2.m G.m
  55. pause % Press any key to graph x = g1(x,y) and y = g2(x,y).
  56.  
  57. clc; clg;
  58. %    Plot x = g1(x,y) and y = g2(x,y) over the rectangle [a,b]x[c,d].
  59.  
  60. a = -2.0;
  61. b =  2.5;
  62. c = -1.0;
  63. d =  1.5;
  64. hx = (b-a)/31;
  65. hy = (d-c)/31;
  66. [X Y] = meshdom(a:hx:b, c:hy:d);...
  67. Z = g1(X,Y) - X;...
  68. W = g2(X,Y) - Y;...
  69. contour(Z,[0 0],a:hx:b,c:hy:d,'-g');...
  70. grid;...
  71. hold on;...
  72. contour(W,[0 0],a:hx:b,c:hy:d,'-g');...
  73. title('Implicit plot of x = g1(x,y) and y = g2(x,y).');...
  74. hold off;...
  75. shg; pause    % Press any key to perform fixed point iteration.
  76.  
  77. clc;
  78. %    Place the starting value in  [p0 q0]
  79.  
  80. % Place the tolerance in  delta
  81.  
  82. %    Place the number of iterations in  max1
  83.  
  84. p0 = 0.0;
  85. q0 = 1.0;
  86. P0 = [p0 q0];
  87. delta = 1e-12;
  88. max1  = 50;
  89.  
  90. [P0,err,P] = fix2dim('G',P0,delta,max1);
  91.  
  92. pause    % Press any key for the list of iterations.
  93.  
  94. clc; clg;
  95. a = -0.275;
  96. b =  0.025;
  97. c =  0.99;
  98. d =  1.002;
  99. hx = (b-a)/20;
  100. hy = (d-c)/20;
  101. [X Y] = meshdom(a:hx:b, c:hy:d);...
  102. Z = g1(X,Y) - X;...
  103. W = g2(X,Y) - Y;...
  104. contour(Z,[0 0],a:hx:b,c:hy:d,'-g');...
  105. hold on;...
  106. contour(W,[0 0],a:hx:b,c:hy:d,'-g');...
  107. X0 = P(:,1);...
  108. Y0 = P(:,2);...
  109. plot(X0,Y0,'or');...
  110. plot([a b],[0 0],'b',[0 0],[c d],'b');...
  111. title('Graphical presentation of the fixed point iteration.');...
  112. grid;...
  113. hold off;...
  114. shg; pause    % Press any key to continue.
  115.  
  116. clc;
  117. Mx1 = 'Computations for the fixed point iteration method.';
  118. Mx2 = '     p(k)               q(k)';
  119. Mx3 = 'The solution is:';
  120. Mx4 = 'The error estimate for P is ';
  121. Mx5 = num2str(err);
  122. Mx6 = [Mx4,'[± ',Mx5,', ± ',Mx5,']'];
  123. clc,echo off,diary output,...
  124. disp(''), disp(Mx1),disp(''), disp(Mx2),disp(P),...
  125. disp('Iteration is converging linearly to the root.'),...
  126. disp(''),disp(Mx3),disp(''),disp('G(P) = P = '),disp(P0),...
  127. disp(''),disp([Mx6]),diary off,echo on
  128.  
  129. pause    % Press any key to perform fixed point iteration.
  130.  
  131. clc;
  132. %    Place the starting value in  [p0 q0]
  133.  
  134. % Place the tolerance in  delta
  135.  
  136. %    Place the number of iterations in  max1
  137.  
  138. p0 = 2.0;
  139. q0 = 0.0;
  140. P0 = [p0 q0];
  141. delta = 1e-12;
  142. max1  = 6;
  143.  
  144. [P0,err,P] = fix2dim('G',P0,delta,max1);
  145.  
  146. pause    % Press any key for the list of iterations.
  147.  
  148. clc; clg;
  149. a = -2;
  150. b =  10;
  151. c = -3;
  152. d =  1;
  153. hx = (b-a)/31;
  154. hy = (d-c)/31;
  155. [X Y] = meshdom(a:hx:b, c:hy:d);...
  156. Z = g1(X,Y) - X;...
  157. W = g2(X,Y) - Y;...
  158. contour(Z,[0 0],a:hx:b,c:hy:d,'-g');...
  159. hold on;...
  160. contour(W,[0 0],a:hx:b,c:hy:d,'-g');...
  161. X0 = P(:,1);...
  162. Y0 = P(:,2);...
  163. plot(X0,Y0,'or');...
  164. plot([a b],[0 0],'b',[0 0],[c d],'b');...
  165. title('Graphical presentation of the fixed point iteration.');...
  166. grid;...
  167. hold off;...
  168. shg; pause    % Press any key to continue.
  169.  
  170. clc,echo off, diary output,...
  171. disp(''), disp(Mx1),disp(''), disp(Mx2),disp(P),...
  172. disp('Iteration is diverging to "infinity".'),...
  173. disp('Note the "scale factor" for the table of computations.'),...
  174. diary off, echo on
  175.